1. Field of the Invention
This method is concerned with automated segmentation of geophysical data for the purpose of providing a model of the shapes and properties of earth layers beneath a region of interest.
2. Discussion of Related Art
Geophysical exploration of the earth's subsurface is commonly done by illuminating the subsurface earth layers using suitable radiation. The radiation may be acoustic, electromagnetic, radioactive or of any other nature now known or unknown. During a survey, a radiation source, positioned at or near the surface of the earth, visits a plurality of survey stations arranged in a grid-like pattern over a region of interest. The source is associated with an array of spaced-apart receivers that are responsive to the radiated energy after that energy has passed through and has been filtered and reflected by the subsurface earth layers. At each survey station, the source is activated to allow the respective receivers to record a snap-shot image of a portion of the subsurface in the light of chosen radiation. Each resulting image is digitally recorded on an archival storage medium in the form of a plurality of time-scale traces, usually one trace per receiver group location. Each trace exhibits a series of transients or waveforms in terms of trace amplitude vs. source-to-receiver travel time. On the multi-trace recordings the respective waveforms of the series may be representative of the stratigraphic sequence of a volume of the earth in the vicinity of each survey station, plus unwanted glare or noise.
It is usual practice to arrange the field procedure to allow redundant overlap of adjacent records. Preliminary routine data processing contemplates filtering out the noise, correcting the data for angularity, spherical spreading and sharpening the wavelet imagery using various well-known stacking and multiple-reflection abatement and migration routines.
From the survey, with particular reference to well-known seismic surveying practices by way of example but not by way of limitation, the processed records are merged to produce a 3-dimensional (3-D) image of a volume of the earth in the selected region. To that end, events (locations of a specified phase of waveform) are tracked continuously across adjacent traces and between records by simple inspection, cross-correlation or other well-known means. Each event, representative of a particular geologic horizon or earth layer, is tracked from trace-to-trace and station-to-station to build the desired model of the subsurface earth layers. The desired model may be purely structural to define the limits of a mineral or hydrocarbon deposit. Concomitantly, using attribute methods, some physical characteristic such as the lithology or porosity of a particular rock layer may be mapped.
The interpreted results of a survey may be shown as a topographic map of one or more specific earth layer interfaces or as a sequence of one or more horizontal frames across a data volume, referred to as time slices, to provide subcrop maps.
In a large survey region, it is prohibitively laborious for a human interpreter to follow a plurality of seismic horizons over thousands of stations. The task is preferably relegated to a computer. Several known computer-aided methods of geological interpretation of geophysical data are now presented by way of review. Such methods fall into two categories based on non-uniform- and uniform-time sampling of seismic attributes.
For the first category, Paulson et al., published a paper on horizon tracking entitled Automatic Seismic Reflection Picking, Geophysics v.33, 1968, p431-40. The method first picks the time positions of either all amplitude maxima or amplitude minima on a seismic trace (i.e. a seismic record at one spatial location). This step creates an array of unequally spaced time positions. Second, it uses maximum normalized cross-correlation criterion to match seismic waveform in a time window centered at picked time on one with that in a possibly time-shifted window centered at picked time on a neighboring trace. The continuous trajectory between reflection times in matched windows on contiguous traces is defined as a seismic horizon. Finally, the method associates with each picked time an attribute called grade, which is the product of the length of a trajectory and mean maximum or minimum amplitude along the trajectory. This method was applied to 2-D seismic data since 3-D data did not exist at the time.
Hildebrand et al. in U.S. Pat. No. 5,251,184 issued Oct. 5, 1993, entitled Method and Apparatus for Finding Horizons in 3-D Data, disclose a process that applies one-trace forward and backward search within specified time window and an amplitude tolerance attribute to track horizons in 3-D data. This method is a modification, of Paulson et al. method, where (a) it iterates between forward and backward trace directions among four orthogonally disposed neighbor traces to verify connectivity among picks, (b) it provides an interactive graphical means of introducing a seed point from where the horizon tracking expands laterally, (c) it keeps a history of trajectory expansion starting from the seed point, (d) it enables the user to edit erroneous trajectories, and (e) it displays a tracked cross-section on a graphical screen. Since the seed input and trajectory editing require human intervention, the horizon tracking with this method sets a practical limit on the number of horizons that may be mapped in a given time.
Having tracked one or more horizons over the region, it is of interest to estimate the dip and strike (or azimuth). Dalley et al., in an article entitled Dip and Azimuth Displays for 3-D Seismic Interpretation, in First Break, v. 7, 1989, p86-95, discuss a method for estimating the dip and strike of an automatically tracked or interpreted horizon. The dip and azimuth are spatial (lateral) attributes associated with the tracked horizon time.
The horizon tracking algorithms have problems. First, the tracked horizon (s) sometimes fail to close a loop even on a locally continuous surface. Problems related to loop closure on continuous surfaces include cycle-skipping due to aliasing in the presence of steep dips and inability to extrapolate the horizon across areas of poor record quality. Second, the Hildebrand method requires a seed point to pick a horizon. That is, a human must tell the computer where to start. In this method, since only one attribute--the maximum or minimum amplitude--is used, noise leads to tracking of pseudo horizons especially in areas of poor record quality. Third, both Hildebrand and Dalley methods provide only sparse information confined to particular interpreted horizons out of an entire volume of possible horizons. Fourth, although all horizon tracking based methods provide structural information, i.e., bulk geometrical shapes of horizons, yet they lack stratigraphic and lithologic information, i.e. detailed subdivisions between boundaries of structural layers and their properties such as porosity or fluids filling the pores. Finally, these methods do not provide any information in regions where horizons cannot be tracked.
For the second category, Taner et al., in a paper entitled Complex Seismic Trace Analysis, Geophysics, v.44, 1979, p1041-63, quantified the character of a seismic reflection through the use of Hilbert Transform to calculate temporal attributes: Instantaneous Envelope, Instantaneous Phase and Instantaneous Frequency. Instantaneous Envelope provides the lithologic component and Instantaneous Phase the structural component of stratigraphy. The first advantage of uniform time sampling based method is that it does not require horizon tracking since attributes are estimated at equal increments of time, unlike the first category methods where attributes are calculated only at unequally spaced horizon picks. The second advantage is that it provides an attribute value even in those regions where a horizon from first category methods does not exist. The disadvantage of Taner's attribute method is that the Instantaneous Phase and Frequency functions become unstable at the minima of the Instantaneous Envelope function. To avoid this disadvantage, Bodine, in a 1984 paper entitled Waveform Analysis with Seismic Attributes, given at the 54.sup.th Annual Meeting of the Society of Exploration Geophysicists, proposed averaging the instantaneous attributes in the neighborhood of the peak Instantaneous Envelope. The averages, termed the Response Envelope, Response Phase and Response Frequency were spread over the entire envelope window that embraces several cycles of the seismic waveform. These parameters were more robust than their Instantaneous counterparts but they lacked the required temporal resolution.
The exploration of many geological features, e.g., faults, meandering river channels, reefs and point bars, is not amenable to either Dalley et al.'s dip and azimuth method or Taner et al.'s temporal attributes method. It requires the mapping of various grades of discontinuity.
In U.S. Pat. No. 5,563,949, issued Oct. 8, 1996, entitled Method of Seismic Signal Processing and Exploration, M. S. Bahorich et al. disclose a method of obtaining a set of discontinuity attribute traces distributed over a 3-D volume of the earth. The 3-D volume of input seismic traces is divided into a plurality of vertically stacked and spaced apart horizontal slices. Each slice is divided into a plurality of cells having at least three seismic traces located therein. The maximum time lagged cross-correlation is measured between the first and second traces lying in a vertical plane to obtain an in-line value and another maximum time lagged cross-correlation between the first and third traces lying in the orthogonal vertical plane to obtain a cross-line value. The geometric mean of the in-line and cross-line values is termed `coherency` associated with the time at the center of the first trace window. Coherency is normalized to lie in the [0,1] interval. The discontinuity attribute is equal to (1-Coherency). This attribute is displayed instead of the seismic time slice on a conventional interpretation workstation. The method works well in areas of gentle dips and low noise, but because it requires an ambient dip and azimuth information it gives false results in areas of steep dips and high noise.
Marfurt et al. at the 67.sup.th Meeting of the SEG, 1997, in a paper entitled Coherency Calculation in the Presence of Structural Dip, Expanded Abstracts p566-7, improved upon Bahorich method by first searching for the best estimate of local dip and azimuth, and then calculating coherence along that dipping plane. The Marfurt technique highlights edges of dipping horizons and enhances features within gaps of horizon tracks in the seismic volume and thus complements the horizon tracking procedure that detects only continuous surfaces. However, the method fails in areas of very steep dips that cause spatial aliasing and provides only one attribute `Coherence`, which reveals only one aspect of underlying stratigraphy.
The basic concept underlying the Marfurt technique is similar to Paulson's method with two exceptions. First, Marfurt's method is a uniform time sampling based method. Here the cross-correlation (or Coherency) value is associated with the time at the center of an equally incrementing time window rather than at a peak or trough as in Paulson's method. Second, it uses normalized cross-correlation itself as an attribute unlike the Paulson method which calculates cross-correlation to connect peaks or troughs on neighbor traces but uses the product of trajectory length and mean amplitude as the attribute. The normalized cross-correlation is a more useful attribute to measure similarity between spatially separated and dipping trace windows.
In U.S. Pat. No. 5,724,309, issued Mar. 3, 1998, W. G. Higgs et al. disclose a method to estimate dip and azimuth from spatial derivatives of Instantaneous Phase. Unlike Dalley et al., it is a uniform time sampling method. In comparison with the first step in Marfurt method, this method directly provides the highly desirable dip and azimuth attributes without a compute-intensive dip search. This method too, however, suffers from false results in areas of high noise or spatial aliasing due to steep dips.
The problems of storage, fast access, and graphical rendering of multiple attribute volumes become more severe with each additional attribute. Oliveros et al., at the 67.sup.th Meeting of the SEG, in a paper entitled Image-Processing Display Techniques applied to Seismic Instantaneous Attributes over the Gorgon gas field, North-West Shelf, Australia, Expanded Abstracts p2064-7, suggest combining 3 attributes by assigning them to govern hue, lightness and saturation in a single color display. This method does not solve the data logistics problem and limits to three the number of attributes that can be visualized at any one time.
To summarize existing known methods for automated interpretation of geophysical data: Horizon tracking requires undesirable extensive human intervention to pick seed-points and edit trajectories. Inconsistent horizon tracks may result due to the use of but a single temporal attribute such as a minimum or maximum amplitude which per se is not uniquely diagnostic of the inter-trace identity of a horizon. Dip and azimuth estimates are restricted to interpreted horizons. Temporal waveform analysis as exemplified by Taner et al. becomes unstable in certain situations and lacks spatial resolution. The Bodine technique of Response Attributes is more stable but it lacks the temporal resolution which is so important in stratigraphic analysis. Coherence and uniformly sampled Dip-Azimuth techniques provide spatial resolution but reveal only one aspect of stratigraphy and fail in the presence of spatial aliasing and noise. Multiple attributes from uniform time sampling suffer from problems of large volume data logistics and visualization. Displays of processed data show either a 3-D volume of a single attribute or a multiple-attribute classification painted in color on a 2-D structural map of a horizon. None of above methods allows a user to take advantage of several attributes and view the joint response in a volume. All methods are fraught with problems of poor loop closure, cycle-skipping in the presence of steep dips and false results in areas of severe noise. In effect, the interpreter of geophysical data must perform many steps, each providing one piece of the structural or stratigraphic puzzle. The interpreter must then synthesize the pieces into a whole picture while remaining aware of the pitfalls.
There is a need for an all-inclusive, cost-efficient, automated process that relieves the user of the need for human-intensive horizon interpretation, editing and mapping. The process should estimate the temporal, spatial and statistical attributes that allow the resolution of the fine structure and stratigraphy of the seismic volume as a whole rather than confining them to a selection of mapped horizons. The process should require modest computer memory and storage resources. And it must offer the user interactive means for speedily and economically exploring the structure and stratigraphy in a "what if" sense having available many combinations of attributes to classify, visualize, and locate hydrocarbon reservoir targets directly in three-dimensional displays.